#include "Multifit.h"

namespace RPI_GSL
{

Multifit::Multifit(int nn, int tt) : n(nn), t(tt)
{
	X = gsl_matrix_alloc (n, t);
	y = gsl_vector_alloc (n);
	work = gsl_multifit_linear_alloc (n, t);
	c = gsl_vector_alloc (t);
	cov = gsl_matrix_alloc (t, t);
	p = new double[t];
}

void Multifit::solve(double *yy, double *xx)
{
	double chisq;
	// fill in the vector and the matrix
	for (int i=0; i<n; ++i)
	{
		gsl_vector_set (y, i, *(yy++));
		gsl_matrix_set (X, i, 0, 1.0);  // x_0 = 1 the constant coefficient
		for (int j=1; j<t; j++)
		{   // set x_1, x_2, ..., x_{t-1}
			gsl_matrix_set (X, i, j, *(xx++));
		}
	}
	gsl_multifit_linear (X, y, c, cov,
                               &chisq, work);
    for (int i=0; i<t; i++){
    	p[i] = gsl_vector_get(c,i);
    }
}

Multifit::~Multifit()
{
	gsl_multifit_linear_free (work);
	gsl_matrix_free (X);
	gsl_vector_free (y);
	gsl_vector_free (c);
	gsl_matrix_free (cov);
}

} // namespace
